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ABSTRACT 


The highest-energy known gamma-ray sources are all located within 0.5 degrees of extremely pow- 
erful pulsars. This raises the question of whether ultra-high-energy (UHE; > 56 TeV) gamma-ray 
emission is a universal feature expected near pulsars with a high spin-down power. Using four years 
of data from the High Altitude Water Cherenkov (HAWC) Gamma-Ray Observatory, we present a 
joint-likelihood analysis of ten extremely powerful pulsars to search for subthreshold UHE gamma-ray 
emission correlated with these locations. We report a significant detection (> 30), indicating that UHE 
gamma-ray emission is a generic feature of powerful pulsars. We discuss the emission mechanisms of 
the gamma rays and the implications of this result. The individual environment, such as the magnetic 
field and particle density in the surrounding area, appears to play a role in the amount of emission. 


1. INTRODUCTION 


Ultra-high-energy (UHE; > 56 TeV) gamma-ray emis- 
sion can be created via hadronic or leptonic processes. 
In the hadronic mechanism, a neutral pion decays into 
two gamma rays. In the leptonic mechanism, a lower- 
energy photon scatters off a relativistic electron via in- 
verse Compton scattering. The electron transfers part of 
its energy to the gamma ray, resulting in a higher-energy 
photon. 

Traditionally, it was thought that UHE gamma-ray 
sources would be hadronic in nature, as leptonic emis- 
sion is suppressed in this energy regime due to the Klein- 
Nishina (KN) effect. However, present-day gamma-ray 
observatories have the sensitivity required to detect lep- 
tonic sources above 56 TeV. See, for example, the de- 
tections of the Crab Nebula as well as several known 
pulsar wind nebulae (PWN) (Abeysekara et al. 2019; 
Amenomori et al. 2019; Abeysekara et al. 2020; Abdalla 
et al. 2019). The astrophysical spectrum of a leptonic 
source typically exhibits significant curvature due to the 
KN effect (Moderski et al. 2005). Conversely, hadronic 
sources follow the spectrum of their parent cosmic-ray 
population, which may or may not include a cutoff or 
curvature. 

The High Altitude Water Cherenkov (HAWC) Obser- 
vatory is a gamma-ray observatory with an wide instan- 
taneous field-of-view (72 steradians) and sensitivity to 
energies between a few hundred GeV and a few hun- 
dred TeV. It is sensitive to sources with declinations 
between -26 and +64 degrees (Smith 2015; Abeysekara 
et al. 2017, 2019). 

The first HAWC catalog of UHE sources (Abeysekara 
et al. (2020); hereafter referred to as the *eHWC" cata- 
log) contains nine sources emitting above 56 TeV, three 
of which continue above 100 TeV. The highest-energy 
sources all exhibit curvature in the spectrum. Addition- 
ally, all nine sources are located within 0.5 degrees of 


pulsars. For eight of the nine sources, the pulsar has an 
extremely high spin-down power (È > 10% erg/s). This 
is much higher than the number of high-E pulsars that 
would be expected to be found near UHE gamma-ray 
sources (Abeysekara et al. 2020). Emission near pulsars 
is expected to be powered by a PWN or TeV halo and 
is therefore dominantly leptonic, even at the high en- 
ergies studied here (Breuhaus et al. 2020; Sudoh et al. 
2021). While young pulsars are expected to have an as- 
sociated supernova remnant (SNR) with accompanying 
hadronic emission, there has only been one detection of 
gamma-ray emission from an SNR to UHE, and leptonic 
emission from this source has not been conclusively ruled 
out (Albert et al. 2020a). SNR theory starts to run into 
technical problems accelerating particles to these ener- 
gies (Gabici 2017). SNR acceleration to UHE energies 
may not be possible (Zeng et al. 2021). 

The proximity of these gamma-ray sources to the 
most powerful pulsars, along with the curvature in their 
spectra, invites the question of whether UHE gamma- 
ray emission is a generic feature expected near these 
sources. This is investigated in this paper through a 
joint-likelihood analysis of pulsars with E > 10% erg/s 
to search for subthreshold sources in the HAWC data. 
While each source is too weak to be individually de- 
tected in HAWC’s standard catalog search, analyzing 
the regions jointly may lead to a general detection for 
this source class. 

'The paper is organized as follows: Section 2 describes 
the details of the joint-likelihood method. Section 3 con- 
tains the results of the analysis. Section 4 discusses im- 
plications of the results. 


2. ANALYSIS METHOD 
2.1. Joint-Likelihood Method 


In this paper, we perform a joint analysis of several 
high-E pulsars. The analysis uses a binned maximum- 
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likelihood method that searches for a gamma-ray excess 
above the background, using a simple model that de- 
scribes the UHE emission from pulsars based on various 
observables such as the distance and pulsar age. This 
is performed using the HAWC Accelerated Likelihood 
(HAL)! plugin to the Multi-Mission Maximum Likeli- 
hood Framework (3ML) (Vianello et al. 2015). We de- 
scribe the background rejection and likelihood method 
in Abeysekara et al. (2019). HAWC’s background 
comes from two main sources: the Galactic diffuse back- 
ground stemming from unresolved gamma-ray sources, 
and cosmic rays that are detected by HAWC and survive 
gamma/hadron separation cuts. Note that neither back- 
ground is very large above 56 TeV. The Galactic diffuse 
background is largely due to Inverse Compton scatter- 
ing, which is not very prominent at these energies. The 
fraction of cosmic rays surviving gamma/hadron sepa- 
ration cuts is also very low: ~0.001 (Abeysekara et al. 
2017). 

Data were collected over 1343 days between June 
2015 and June 2019. The data are binned in quarter- 
decade bins of reconstructed gamma-ray energy using 
the “ground parameter” energy estimator, one of two 
energy estimators currently used by HAWC (Abeysekara 
et al. 2019). 

We use the last three quarter-decade energy bins from 
Abeysekara et al. (2019), restricting the analysis to re- 
constructed energies between 56 and 316 TeV. These are 
the highest energies probed by HAWC. At these ener- 
gies, there is very little diffuse emission and less source 
confusion than at lower energies. 

'The spectral model is a power law: 


Ear (E). (1) 


where A; accounts for the model-dependent relative flux 
of each source (see Section 2.3), K is the normalization, 
and Ep is the pivot energy, which is fixed at the cen- 
ter of each energy bin. The three values used for Ep 
in this analysis are 74.99 TeV, 133.35 TeV, and 237.14 
TeV. Setting the pivot energy at the center of each bin 
makes the analysis relatively insensitive to the choice of 
a, which is fixed at 2.5. 

We fit each source (see Table 1) individually in each 
energy bin, placing a step function at the boundaries of 
the bin to ensure that events that are mis-reconstructed 
in energy are excluded. We assume the sources are spa- 
tially extended with a Gaussian morphology. The width 
is fixed at 0=0.34%. The values of a and o are approxi- 


l https://github.com/threeML/hawc.hal 


mately the average values for the highest-energy gamma- 
ray sources from the eHWC catalog (Abeysekara et al. 
2020). 

'To obtain a joint-likelihood result, the log-likelihood 
profiles for each individual source, in each energy bin, 
are added and the value of K (hereafter called K) that 
optimizes this log-likelihood profile is found. The total 
flux normalization, &, for the ten sources combined is 
then just: 


Nsources 


wee Ye BA (2) 
i=1 
Using this flux normalization, the total differential 
flux from all ten sources is: 


dN EN? 
ew. 3 
dE ^ (5.) (5) 
We calculate a test statistic (TS) to show how signif- 
icant the value of k is: 


TS = zi ZBC), (4) 
Lp 

where Ls+8 is the maximum likelihood for the signal- 

plus-background hypothesis and Lg is the likelihood for 

the background-only hypothesis. 

After an overall best-fit value of k is determined, it is 
used as an input to a Markov chain Monte Carlo and 
a distribution of the value is determined. If the TS in 
a given energy bin is significant (TS > 4), a Bayesian 
credible interval (68% containment) is shown, obtained 
from the estimated distribution of the parameter &. The 
model only has one degree of freedom, so using Wilks 
Theorem (Wilks 1938), a TS value of 4 corresponds to 
2c. Otherwise, 95% credible interval quasi-differential 
upper limit is plotted. A uniform prior is used in the 
Bayesian analysis. 

In the bins where an upper limit is determined, the 
range of expected upper limits are obtained from simu- 
lations of Poisson-fluctuated background. 


2.2. Source Selection and Dataset 


For this analysis, we define the source list by select- 
ing all pulsars from the ATNF pulsar database, version 
1.62? (Manchester et al. 2005) in the inner Galactic 
plane that are within HAWC’s field-of-view (|b| « 1°, 
5° < | < 90°) and have E > 10% erg/s. There are 24 
pulsars that meet this criteria (see Table 1). Pulsars 
with a high E are centered around the Galactic plane, 
so the choice to concentrate on |b| « 1? only removes 


2 https://www.atnf.csiro.au/research/pulsar/psrcat/ 
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Figure 1. VTS for É > 56 TeV emission from HAWC, for the field-of-view used in this analysis. We assume a point source 
morphology. The locations of the ten sub-threshold pulsars used in this analysis are labeled. 


3 additional pulsars from the analysis. We use this list 
of 24 pulsars to determine which models of emission are 
reasonable (see Section 2.3). 

After making this determination, this list is then 
downselected to search for sub-threshold gamma-ray 
sources. First, we remove all pulsars that are lo- 
cated within a degree of sources from the eHWC cata- 
log (Abeysekara et al. 2020). This removes pulsars that 
already have significant UHE emission detected in their 
vicinity. Since the gamma-ray emission is modeled as 
extended in nature, this also removes pulsars whose as- 
sociated emission may overlap with those known sources, 
which would require more detailed modeling. 

We remove three additional pulsars from the source 
selection. PSR J1813-1749 is removed because HAWC 
has added more data since the publication of Abey- 
sekara et al. (2020) and has detected a new UHE source 
(eHWC-J1813-176) ~0.2° away from the pulsar. PSR 
J1913+1011 and PSR J1930+1852 have been removed 
because the known TeV emission in those regions is 
likely from a SNR, not a PWN (Abdalla et al. 2018a,b). 
Since the majority of the emission at these energies is 
expected to be from a PWN or TeV halo, this is done to 
prevent the introduction of a different, predominantly 
hadronic source class into the analysis. 

The final list is composed of ten pulsars that are can- 
didates for sub-threshold gamma-ray emission. Figure 
1 shows HAWC’s > 56 TeV map with these sources la- 
beled. 


2.3. Models 


The parameter A; in Equation 1 describes the rela- 
tive contribution each pulsar receives in the analysis. 
This parameter can be used to test different models 
of gamma-ray emission near pulsars. For the models 
that rely on pulsar parameters, the relevant quanti- 
ties are taken from the ATNF pulsar catalog, version 
1.62 (Manchester et al. 2005). We consider a variety of 
different models. In all descriptions, “emission” refers 
solely to gamma-ray emission above 56 TeV. 


e No model: Here, A; for all sources is set equal to 
1. All sources are treated equally and the emission 


is expected to be uncorrelated with pulsar param- 
eters such as distance. 


e 1/d?: In this model, A; is set to 1/d?, where d 
is the distance to the pulsar. This model assumes 
that closer sources will produce observable emis- 
sion. 


e E/d?: Here, the 1/d? model discussed above is 
multiplied by the spin-down power of the pul- 
sar. Therefore, closer, more-energetic pulsars have 
more gamma-ray emission. 


e Inverse age: In this model, A; is the inverse of 
the spin-down age. This is defined as P/(2P), 
where P and P are the period and the time deriva- 
tive of the period, respectively. In this model, 
younger sources are more likely to have detectable 
emission. 


e Flux at 7 TeV: Here, A; is the HAWC flux 
at 7 TeV. This model assumes that sources that 
are bright at multi-TeV energies should also give 
off detectable emission above 56 TeV. The 0.5? 
extended source map from the third catalog of 
HAWC sources (3HWC) (Albert et al. 2020b) is 
used to extract these values. A; is computed by 
averaging the flux from all pixels within a 0.5 de- 
gree radius of the pulsar. 


3. RESULTS 


We first run the joint-likelihood analysis (as described 
in Section 2) with the full list of 24 pulsars to determine 
which models of emission are reasonable. Table 2 gives 
the total TS for each scenario. For each model, we can 
also calculate the expected gamma-ray flux above 56 
TeV from an arbitrary pulsar using Equation 1. The 
values of K, derived from Equation 2 are also given in 
Table 2. In all cases, the TS is much higher than 50. 
Two scenarios, the flux at 7 TeV and 1/d?, perform 
better than the “no model case". 'The model based on 
the inverse age of the pulsar performs significantly worse 
than the others, but the significance, VTS is still well 
above the 5c level. 
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PSR name | RA (°) | Dec (^) | Age (5) (kyr) | E (x10% erg/s) | Distance (kpc) | P (s) P (ss ^1) Subthreshold? 
B1800-21 270.96 | -21.62 15.8 2.2 4.40 0.134 | 1.35 x10? Y 
J1809-1917 || 272.43 | -19.29 51.3 1.8 3.27 0.083 | 2.55 x107** 
J1811-1925 272.87 | -19.42 23.3 6.4 5.00 0.065 | 4.4 x10714 
J1813-1749 || 273.40 | -17.83 5.6 56 4.70 0.045 | 1.27 x10713 
J1826-1256 276.54 | -12.94 14.4 3.6 1.55 0.110 | 1.21 x10? 

B1823-13 276.55 | -13.58 21.4 2.8 3.61 0.101 | 7.53 x107** 
J1828-1101 277.08 | -11.03 77.1 1.6 4.77 0.072 | 1.48 x10 * "4 
J1831-0952 || 277.89 -9.87 128 1.1 3.68 0.067 | 8.32 x107* Y 
J1833-1034 || 278.39 | -10.57 4.85 34 4.10 0.062 | 2.02 x10? "4 
J1837-0604 || 279.43 -6.08 33.8 2.0 4.77 0.096 | 45.1 x107* 
J1838-0537 | 279.73 -5.62 4.89 6.0 2.0* 0.146 | 4.72 x10 
J1838-0655 || 279.51 -6.93 22.7 5.5 6.60 0.070 | 4.93 x107** Y 
J1844-0346 281.14 -3.78 11.6 4.2 2.40” 0.113 | 1.55 x107*8 
J1846-0258 || 281.60 -2.98 0.728 8.1 5.8 0.327 | 7.11 x10? 
J1849-0001 282.23 -0.02 42.9 9.8 7.0° 0.039 | 1.42 x10 * 
J18564-0245 | 284.21 2.76 20.6 4.6 6.32 0.081 | 6.21 x107** "4 
J19074-0602 | 286.98 6.04 19.5 2.8 2.37 0.107 | 8.68 x107** 
J1913+1011 | 288.33 10.19 169 2.9 4.61 0.036 | 3.37 x107" 
J19284-1746 | 292.18 17.77 82.6 1.6 4.34 0.069 | 1.32 x107** Y 
J19304-1852 | 292.63 18.87 2.89 12 7.00 0.137 | 7.51 x10713 
J1935+2025 | 293.92 20.43 20.9 4.7 4.60 0.080 | 6.08 x107** Y 
B1937+21 294.91 21.58 2.35e5 11 3.50 0.002 | 1.05x10- !? Y 
J20214-3651 || 305.27 36.85 17.2 3.4 1.8 0.104 | 9.57 x107** 
J20224-3842 | 305.59 38.70 8.94 30 10.00 0.049 | 8.61 x107** Y 


?Pseudo-distance from Pletsch et al. (2012) 
PPseudo-distance derived from Eq. 3 of Wu et al. (2018) 
*Distance estimate from Gotthelf et al. (2011) 


Table 1. Information on the pulsars. All information comes from the ATNF database, version 1.62 (Manchester et al. 2005) 
except for some distance estimates which are not included in the pulsar database (see footnotes). Here, E is the spin-down 
energy loss rate, P is the barycentric period, and P is the time derivative of the period. The checkmark in the last column 
denotes the ten pulsars included in the sub-threshold analysis. 


Model K (56 < E < 100 TeV) | K (100 < E < 177 TeV) | K (177 < E < 316 TeV) Units for A; Total TS 
No model 2.55 x 107* 4.65 x 107" 9.46 x 10718 dimensionless 633 
1/d? 244 x 10718 4.50 x 10-16 8.12 x 10 kpc ? 734 
E/@ 4.06 x 107% 6.88 x 10 + 1.04 x 10 * ergm ? s! 568 
Inverse age 6.36 x 10718 1.10 x 10718 1.65 x 107 4 yr! 152 
Flux at 7 TeV 5.42 x 10? 9.11 x 1074 1.90 x 1074 TeV! cm”? s^! 886 


Table 2. The proportionality constants needed to calculate the expected gamma-ray flux near a given pulsar. The dimensions 
for A¿K are energy! distance ! time !. With the value of K known and A; in the units given by the last column, the expected 
gamma-ray flux can easily be calculated using Equation 1. K is reported at the pivot energy in each bin. Proportionality 
constants are desrived from the full sample of 24 pulsars. 
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Figure 2. The total combined flux from the ten sub- 


threshold candidates for the scenario where the pulsars have 
a relative contribution defined by 1/d?, where d is the dis- 
tance between the pulsar and the Earth. For the first two 
energy bins, the TS > 4 so Bayesian credible intervals (68% 
containment) are plotted. In the last energy bin, there is no 
significant detection so a 95% upper limit is plotted. The 
dark and light grey bands are 68% and 9096 containment for 
expected upper limits from Poisson fluctuated background. 
'The dotted lines are systematic uncertainties on the central 
value (i.e., the solid red line may be as low or as high as the 
dotted line once systematic uncertainties are included, see 
section 3.1) 


Satisfied that the chosen models are adequate, we then 
run the joint-likelihood analysis using the downselected 
list of ten pulsars (as described in Section 2.2) to search 
for sub-threshold leptonic emission from pulsars. Re- 
sults are given in Table 3. Each of the models give a 
total TS value above 9, corresponding to a detection of 
more than 3c. For some models, the TS is much higher. 
The 1/d? model gives the highest TS: 41.3 (6.40). The 
same two models as before, 1/d? and the gamma-ray flux 
at 7 TeV, perform better than the ^no model” case, but 
in this case 1/d? performs the best, with a TS slightly 
higher than the “gamma-ray flux at 7 TeV” scenario. 

Figure 2 shows the total flux from all ten sub-threshold 
candidates for this best-case scenario. The figures for 
the other models can be seen in the Appendix. Table 4, 
also located in the Appendix, gives the total flux nor- 
malizations for each model. 

'The Appendix also includes a discussion of how often 
TS this high would be expected from stacking random 
points on the sky. 


3.1. Systematic Uncertainties 


Systematic uncertainties are broken down into two 
categories: detector systematics and modeling system- 
atics. Detector systematics, described in Abeysekara 
et al. (2019), stem from mis-modeling of detector quan- 
tities such as the photomultiplier tube threshold and 
charge in simulated Monte Carlo events. Each system- 


atic is treated independently; the results are added in 
quadrature to get a total uncertainty. This process is 
repeated in each energy bin. Depending on the energy 
bin and model assumed, the detector systematic ranges 
from 1096 to 2596. 

Modeling systematic uncertainties investigate how the 
analysis would change if the emission near the sub- 
threshold pulsars is different from what is assumed in the 
main analysis. Several modeling systematics are consid- 
ered: 


e Spectral indices in the power law (a in Equation 
1) of 2.0 and 3.0 


e Replacing the power-law spectral model (Equation 
1) with a power-law with an exponential cutoff: 


dN, EN —E/Ecut 
AK (57) e . (5) 


Here, a is fixed at 2.5 and E, is fixed at 60 TeV, 
which are average values from Abeysekara et al. 
(2020) 


e Changing the Gaussian width to 0.23? and 0.45”. 
These are + 1 standard deviation from the aver- 
age extension of the sources from HAWC's eHWC 
catalog (Abeysekara et al. 2020). 


'The modeling systematics are larger than the detec- 
tor systematics, driven predominantly by the assumed 
source size. In the last energy bin, assuming a power 
law with an exponential cutoff instead of a power law is 
also a dominant effect. 

Depending on energy and model, the modeling sys- 
tematic ranges from 13% to 34%. The sum of the de- 
tector and modeling systematics, added in quadrature, 
are denoted with a dotted black line in Figure 2 and in 
Figures 3 through 5 in the Appendix. 

As described in Albert et al. (2020b), the absolute 
pointing uncertainty of HAWC is declination-dependent 
and may be larger than previously thought; perhaps as 
large as 0.3 at the edges of the field-of-view. This means 
that the TS values presented in this paper may be un- 
derestimated; the size of this effect is ~10%. The flux 
may also be underestimated. The flux underestimate is 
a function of energy and ranges from 9% in the 56-100 
TeV bin, decreasing to a negligible amount (0.6%) in the 
177-316 TeV bin. 


4. DISCUSSION 


Regardless of which model is used, it is clear that the 
areas around high-E pulsars show hints of emitting at 
ultra-high energies (> 56 TeV). Interestingly, some mod- 
els based on the pulsar parameters, such as E or inverse 
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Model TS (56 « E « 100 TeV) | TS (100 < E « 177 TeV) | TS (177 « E « 316 TeV) | Total TS 
No model 27.9 8.33 1.59 37.9 
1/d? 31.9 9.08 1.29 41.3 
E/d? 9.58 5.24 0.00 14.8 
Inverse age 9.19 3.79 0.03 13.0 
Flux at 7 TeV 26.3 9.61 3.62 39.6 


Table 3. The test statistic for the joint-likelihood analysis for each model, using the ten sub-threshold sources. Note that « is 
fit individually in each energy bin so the last column is not a joint TS. 


age, have much lower TS values than the “no model” 
case, implying that they do not describe the emission 
well. Two models that perform better than the “no 
model” case are the gamma-ray flux at 7 TeV and 1/d?. 

It is unclear at this time why some pulsar parameters 
do not seem to be good predictors of emission. Some of 
these parameters have fairly large uncertainties, which 
could be a contributing factor. For example, the char- 
acteristic age of the pulsar can differ from the true age. 

However, the uncertainties on the pulsar distances are 
relatively small. For the ATNF pulsar database, 95% 
of pulsars will have a relative error of less than a factor 
of 0.9 in their distance estimate (10% uncertainty) (Yao 
et al. 2017). 

Alternatively, the individual environment that each 
pulsar is in could play a large factor in the amount of 
UHE emission. While the pulsar itself is the particle ac- 
celerator, diffusion of electrons and positrons away from 
the pulsar is strongly dependent on quantities such as 
the density and magnetic field of the environment. Note 
that this may be only true for the high-E pulsars stud- 
ied here, which are relatively young; this means the mag- 
netic fields are likely to be affected by the pulsar age and 
SNR interactions. For weaker, older pulsars, the mag- 
netic field density and environment are instead likely 
dominated by the ISM. 

Also, the emission mechanisms are still uncertain. 
While it is commonly assumed that the bulk of emis- 
sion from the vicinity of a pulsar is from a PWN and 
therefore predominantly leptonic, a hadronic contribu- 
tion cannot be a priori discarded. While emission from 
associated SNRs are unlikely, some have raised the pos- 
sibility of more exotic hadronic emission mechanisms in 
or near PWN (Amato et al. 2003; Di Palma et al. 2017). 
Several of the HAWC sources known to emit above 56 
TeV, most notably eHWC J19084-063 and eHWC J1825- 
134, have molecular clouds nearby (Duvidovich et al. 
2020; Voisin et al. 2016). These molecular clouds may 
be serving as a target for a portion of the gamma-ray 
production. 


Multi-wavelength and multi-messenger campaigns can 
help disentangle emission mechanisms. A neutrino de- 
tection coincident with one of these pulsars would be 
a smoking gun for hadronic emission mechanisms in or 
near PWN. However, a recent stacked analysis looking 
for neutrino emission from PWN by IceCube did not 
yield a detection (Aartsen et al. 2020). 

Electromagnetic multi-wavelength studies could also 
be helpful. A leptonic source emitting above 56 TeV 
will have a different signature at lower energies than a 
hadronic one. For example, 100 TeV gamma rays imply 
a synchrotron peak in the keV regime, assuming a 3 
uG field (Hinton & Hofmann 2009), with the emission 
extending up to the MeV energy range. If the emission 
is instead predominantly hadronic, there will be no such 
peak at these energies. Proposed experiments such as 
AMEGO (McEnery et al. 2019) will be important in 
distinguishing emission mechanisms. 


5. CONCLUSIONS 


In this study, we have searched for UHE gamma-ray 
emission in the vicinity of pulsars with an E > 10% 
erg/s. We find, with high significance (> 3c) , that UHE 
gamma-ray emission is a generic feature in the vicinity 
of this class of pulsars. 1/d? is the model that gives the 
highest TS; a source that is closer to us is more likely to 
have observed UHE emission. Other pulsar parameters 
do not seem to be good predictors of emission. This 
implies that the environments the pulsars are located in 
may play a role in the amount of emission. 

'The TS values obtained are higher than would be ex- 
pected from combining random points in the Galactic 
plane. There is the possibility that all known gamma- 
ray sources emit above this energy threshold, albeit at 
an extremely low flux. 

Multimessenger and multiwavelength studies are 
needed to disentagle the origin of the UHE emission 
from these pulsars. 
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APPENDIX 


.1. Additional joint-likelihood results 


Figures 3 through 5 are analogous to Figure 2 in the main text, but for the other four models that have been 
investigated (described in Section 2.3). All figures contain the combined flux for the joint-likelihood analysis of the 
ten sub-threshold pulsars. Table 4 contains the total combined flux normalizations for each of the models that have 
been considered. 


6x 107% 


4x107? 


3x10? 


Flux (TeV / cm^2/ s) 


2x10? 


6x10! 10? 2x10? 3x10? 
Energy (TeV) 


Figure 3. Identical to Figure 2 from the main text, but for the ^no model" case. 


4x1077? 


3x 10-4 


2x10-¥ 


Flux (TeV / cm^2/ s) 


10-22 


6x 10-2 


6x10! 10? 2x10? 3x10? 
Energy (TeV) 


Figure 4. Identical to Figure 2 from the main text, but for the E/d? model. 


4x107”? 


3x1077 


2x10? 


Flux (TeV / cm^2/ s) 


6x10! 10? 2x10? 3x10? 
Energy (TeV) 


Figure 5. Identical to Figure 2 from the main text, but for the inverse age model. 
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6x107? 


4x107”? 


3x10? 


Flux (TeV / cm^2 / s) 


2x102 


6x10! 10? 2x10? 3x10? 
Energy (TeV) 


Figure 6. Identical to Figure 2 from the main text, but for the model defined by the gamma-ray flux at 7 TeV. 


Model k (56 < E < 100 TeV) | « (100 < E < 177 TeV) | & (177 < E < 316 TeV) 
No model Sart x 1071$ lac ggs 7.56 x 1071" 
1/d? CU x 10716 154795 x 10716 7.07 x 1071" 
E/@ 3.45+1:27 x 10716 7.414549 x 10717 3.39 x 107 
Inverse age 3.991159 x 10716 1.59 x 10-19 4.48 x 1077 
Flux at 7 TeV 7.807758 x 10718 1.5649: x 10718 7.98 x 1071" 


Table 4. The total combined flux normalization for the sub-threshold joint-likelihood analysis for each of the models (« from 
Equation 3). The units are TeV^! cm”? s^!. The flux normalization is reported at the pivot energy, which is the center of each 
bin. Values without uncertainties are upper limits, otherwise the values correspond to the 6896 containment for the Bayesian 
credible interval. Uncertainties are statistical only. 


.2. Testing of Random Backgrounds 


We investigate how often randomly chosen non-source locations give TS values as high as in the study presented in 
the main text. Sets of ten randomly chosen points from the analysis region (4? < | < 90°; |b| < 1°) are run through the 
joint-likelihood analysis. Points that are within a degree of either a known high-energy source or any of the selected 
A'TNF pulsars are excluded. Because the gamma-ray emission is assumed to be spatially extended, this is necessary 
to avoid contributions from the known sources and the pulsars of interest. 

This analysis is performed for 2,877 sets of ten random source locations. Only five of these randomly chosen sets 
have a TS greater than 37.9. This is the TS from the sub-threshold joint-likelihood “no model” case and is used as 
the comparison because the other models require known pulsar information. This means that the results that we have 
obtained are significant at the 30 level. 

We also investigate sets of ten randomly chosen sources from HAWC’s third catalog (3HWC) (Albert et al. 2020b) 
to see if any of HAWC’s previously detected TeV sources emit at UHE. Once again, sources that are known to emit 
above 56 TeV and sources within a degree of the ATNF pulsars used in the nominal analysis are removed. Note that 
the majority of HAWC sources are leptonic in origin (Linden et al. 2017). 

Due to the relatively small number of remaining 3HWC sources after this downselection (48 of the original 65 3HWC 
sources), it is not possible to run thousands of sets of ten sources without repeating a large subset of the sources. We 
instead run 100 sets of ten randomly chosen 3HWC sources. 

None of these 100 trials have a TS above the value from the sub-threshold joint-likelihood unmodeled case. The 
highest TS from this study is 26.9 and the mean is 12.2 (standard deviation 5.5). This implies that known gamma-ray 
sources that are not located near high-E pulsars are unlikely to emit at high energies, or if they do, their fluxes are 
very low and combining ten sources is not enough for a detection. 

We also explore combining 20 3HWC sources at a time, instead of the ten that were combined in the preceding 
paragraph. The entire TS distribution is shifted to higher values. The highest TS is 45.9 (higher than in the “no 
model” case); the mean is 24.6 and the standard deviation is 7.2. The higher TS values when more 3HWC sources are 
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combined implies that UHE emission may be a generic feature of known gamma-ray sources, but with a very low flux. 
This should be investigated, but may be hard to do with current-generation experiments. Proposed experiments such 
as SWGO (Huentemeyer et al. 2019), will be important here. 


